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1 .0 INTRODUCTION . 

This is the second report on data processing techniques useful for 

, — 

infrared astronomy data analysis systems. As with the first report 
(NASA CR-151943), the investigation is restricted to consideration 
of data from space-based telescope systems operating as survey 
instruments. In this report the theoretical background for specific 
point-source detection schemes is completed* and the development of 
specific algorithms and software for the broad range of requirements 
is begun. 

Section 2 develops the detail detection tests and processing 
requirements for point-source surveys and evaluates the performance 
measurement processes. The details of peak detection decisions 
and correlation detection are covered for the case of general 
bandlimited white gaussian noise. For non-white noise, a modified 
correlation test and a matched filter test are presented. A 
technique for resampling the data which is equivalent to a matched 
filtering approach is discussed which automatically decorrelates 
the noise. Implementation of this Karhunen-Loeve filtering is 
necessarily complicated, but for some kinds of noise an 
acceptable approach. 

Section 3 then reviews a basic processing task to indicate where 
computation is needed outside of the normal data stream. While the 
processing used in the primary data reduction task is important, 
the actual software depends heavily on the specific mission hardware 
and is best approached anew for each task using the theories of 
Section 2 of this report and of the previous report, and of several 
cited authors. For the general signal processing task, a routine 
. for designing digital filters is given based on the theory of 
Section 2.5. The calibration of detector- filter systems is the 
most complicated of the tasks off the main processing line; a routine 
which provides this calibration for blackbody or other input spectra. 
Finally, the preliminary processing routine for a previous survey 
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program is presented and briefly discussed to indicate how much 
processing can be done in a single pass of the data. 

f 

The Appendix in Section 4 presents an interesting game which can 
develop a .fuller appreciation and understanding of the complexities 
of data analysis. 
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2.0 TECHNICAL ASPECTS 


This section completes the task begun in the first report of 
reviewing the theoretical basis for the design of .point-source 
survey data analysis softvvare. The detection techniques for 
single-channel signal and noise processing are reviewed. The 
schemes reviewed include peak detection, optimal filtering, 
correlation, and Karhunen-Loeve filtering. The details of . 
digital filtering, vvhich is applicable to many aspects of data 
processing, are reyiewed-in the final section. 


2.1 Detection of Signals in Noise 

In most coFTimuni cation systems the errors {false detections and 
missed signals) are assumed to be of equal importance and with 
known probabilities. In more general detection problems, however, 
the a priori probabilities and costs of those errors are difficult 
to determine. The Neyman- Pearson test was first applied in such a 
case to radar detection with a peak measurement technique. The 
criterion can also be applied to more sophisticated detection 
methods, and in all cases, will give the highest probability of 
detection at a chosen false-alarm rate. The type of technique used 
depends on the amount of information available about the expected 
signal; generally, more- information used will result in a higher 
detection probability at the chosen false-alarm rate. The likelihood 
ratio is the test used where the hypothesis is chosen if: 


X 


p(s) 


>1 n 


2 . 


and the counter-hypothesis (nosignal) is chosen otherwise. Here 
p(s) is the probability density function of the data with a signal 
present and p(n) is the p.d.f, of the noisealone, and n is the • 
decision level chosen to satisfy the false-alairm constraint. 


Consider the case of a signal in white noise, such that the signal 
has a normalized mean value of one. The probability functions are; 
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2 . 1-2 


p(s)’=_L e-(y-l)V4 and p(n) = J- 

/47r /4 tt 

Then the likelihood ratio test is: 

X(y) = 21 Ap 2.1-3 

To determine the threshold a^, the false-alarm probability is 
found from: 

P(f.a.) = f-L dy . 2.1-4 

Y 

If we want a false-alarm -rate of 10% or less, then y = 1 .8, and v/e 
choose the hypothesis if y >. 1.8. 


The probability of detection for a single test observation is; 


P(det.)= f — dy = 0.285 

■’ /4ir 
Y 


2.1-5 


In terms of the likelihood ratio, note that A(y) = Aq =1.9 and we 
make a detection whenever A(y) ^ 1.9. To improve this rather 
mediocre performance, several measurements may be tested. With 
the same false-alarm rate, we choose the decision level differently. 
If we take n independent samples, the signal-present probability 
distribution has unity mean and a variance of 0 ^, and: 


?s(yi > ^ 2 .* ‘ ’ -^n^ " 


1 

exp [- 

(y^-D^ 

v^T7 a 

2a^ 

1 

exp [- 


cr 

2a^ 


] X . . 


2 . 1-6 


Similarly, the noise-only probability distribution is: 





2.1-7 
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Taking the logarithm of the likelihood ratio, the decision test is: ' 


1 ^ 

— 5] y,. > x' 
n ^^1 •'1 — 0 


2.1-8 


where 


,l_1,a^ Inl 


Nov/ the probability of each kind of error is different, and vie 
choose by evaluating P-j (f.a.) + {missed signal), where 

( , 

P,(f..a.) = I" >„(y) d7= I* ( 5 ^ 2 )'''^ 




and 


P2(m.s.) 


f° P^l.7> dy = e-n(z-l)V 2 o 2 


2.1-9 


It now becomes clear that improving the performance, of simple peak- 
detection schemes becomes a complicated task even using very little 
information about the signal. If we use more of the information 
available, and some of the knowledge about the nature of the noise, 
a more successful detection scheme can be derived. When multiple 
detections- are made on a single source, the above can be used to 
evaluate the detection probability. 

2.2 Correlation Detection 

* C* 

Rather than make a detection test based on only the peaks of the data 
stream as in the previous example, consider how vie migh.t deal with 
detecting a signal that we know. Let r^,, k=l , . . . , m be 
the sequential' data samples. Assuming additive noise, 

=' A} + 2.2-1 
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■Hh 

where is the value of our expected signal, and n^, is the 
noise sample. We may now derive a likelihood ratio test which 
uses this information. 


First, assume that the noise is bandlimited white noise with pov/er 
spectral density S{(i3)' = N^^/Z-for |{ol < and zero otherwise.- The 
noise autocorrelation function- then is given by: ■ 


R(t) 


^00 sin (Qt) 

2tt 


2 . 2-2 


This has its first zero at t = tt/0 so that i.f the received signal 
is sampled at intervals At = ir/Q the samples will be uncorrelated, 
and being gaussian they then will be statistically independent. 

The probability density functions of the two cases will be: 


Ps(i-) 


1 m/2 m 

(^1 .,[- 1 , 


and 


2.2-3 


, m /2 m ( r ,^)2 

P„(r) = (i) exp [- y, ^ ] 


k-1 

and the logarithm of the likelihood* ratio test results in the 
decision test: 


m 

I 

k=l 


’^K^K 


In 




m 

I 

k=l 


S 2 


2.2-4 


Now the left-hand side of 2.2-4 is just the normalized cross- 
correlation coefficient of the signal v</ith its expected template. 
Furthermore, the variance of the noise is just the noise auto- 
covariance function at zero frequency. 
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Since one of our two signals is zero {noise only), we may. define 
the average signal energy E and the time cross-correlation coefficient 
P by: 



m 

k=l 



2 , 2-6 


and 


p = 0 


By extending equations 2.2-3 to infinite, bandwidth the probability 
density functions for the signal case and the noise-only case can be 
derived as: 


1 

2N^E ^ 


1 1/2 /p n\2 


2.2-7 


Since the false-alarm rate and the missed sources probabilities are 
equal when the samples are uncorrelated, the error rate is: 


P 


e 


where 



2 . 2 - 8 . 


Y = 

and thus we can determine an error rate based on the signal-to-noise 
ratio, independent of the shape of the signal. Figure 1 shows the 
error rate as a function of the signal-to-noise power ratio. Note 
that as long- as the noise samples are uncorrelated, the error rate 
is also independent of the number of samples in the correlation 
sum. This apparently unreasonable- result is directly related to 
the assumption of statistically independent samples. For bandlimited 
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white noise there must be m = S2t/ir = constant independent samples in the- 
interval 0 to T. Oversampling the signal may actually result in degraded 
.performance, as will be discussed in section 2.4. If we choose the 
value of Y in 2.2-8 to achieve a desired error rate, then the probability 
of detection is 




- 1/2 
(2ir). 

Y-(2E/Nq)^/^ 


exp [' 


z2/2] dz 


2.2-9 


which is shown in Figure 2 as a function of signal-to-noise ratio and 
error rate. 


If we have chosen the normalized signal template properly, oun detection 
test simultaneously makes a best estimate of the signal amplitude. 

If the signal model is written as a function of a constant amplitude 
factor A, then the maximum likelihood estimate of that amplitude is 
the solution of: 

j [r, - (A)] -2^ = 0 . 2.2-10 

•or, writing s = A s^ , we want the solution of: 


m 

I {r. - A S.)S. = 0 .2.2-11 

i=l ’ ^ ^ 


•That solution is 


m 


A 


m 


i=l ^ 


2.2-12 


_ 2 

and now: if S. was normalized such that > S. = 1,- and we re-arrange 
1 1 . 

the terms in 2.2-4, we have the detection test -and amplitude estimate 
simultaneously: 



2.2-13 
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Now-' it is clear how the correlatioi; test is a better detector than 
a peak test. The correlation test, takes an average of the signals 
weighted by the expected response ^s a best esimate of the amplitude. 
Because it is using n samples of the signal, the improvement in 
error rate can be as much as vfT. The uncertainty in the estimate 
is determined from the noise autocovariance function, as in 
section 2.3 as: 



Rj^(z-t) dzdt 


2.2 


If we have N multiple pulses available from a single source, the 
decision test 2.2-4 can be modified to: 


N m 

I I 

i=l k=l 


^k ^"n 


In X + 


N 

I 

i=l 


m 

I 

k-1 


The false-alarm rate given by 2.2-8 is not changed, since v/e are 
designing our test for a chosen error performance. However, the 
detection probability improves; the new detection rate is given 
by 2.2-9 by replacing E with E': 


N 

E’ = I E. 
i=l ^ 


And since all signal energies for a given source are equal, there 
will be a 3 dB increase in the equivalent performance for each doubling 
of the number of signals. 


As a final note to this discussion, the signal-to-noise ratio used 
here is the more useful signal power- to-noise power ratio, not the 
typical peak-to-rms value which has little physical meaning. “ 



2.3 Matched Filters and Non-White Noise 

The cross-correlation term on the left-hand side of equation 2. 2-4 can 
be replaced by the equivalent matched filter. - If the filter's transfer 
function is h(t^), then the output of the filter is 

m 

■ ,1 

and by inspection the filter -output matches the correlator output 
if 


That is, the .matched filter is the time-reversed image of the signal 
expected. I-t is important to note, however, that because of the time 
reversal, the matched filter and the correlator output are equal 
only at time T, where the entire signal’ train (in samples) is within 
the bounds of the filter or of the correlator. 

The matched filter representation is well suited to the case of non- 
white noise. We will show that the optimal detector for non-white 
noise replaces the left-hand side of 2.2-4 vnth a filter which is the 
product of the white-noise matched’ filter and a pre-whitening filter 
described in terms of the autocovariance function of the noise. To 
avoid confusing subscripts, we shall write the filter transformations 
in terms of time integrals which are the generalized extensions of 
the summations in section 2.2. The output of the filter at time T 
is: 

' e(T) = I h(r) r (T-t) dx ■ 2.3-2 

0 

The signal and noise components are easily. identified as 
S(T) = |\(t) S(T-t) dx 
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h(x) n (T-x) dx 


The noise power then can be v/ritten in terms of the autocovariance 
function as 


a 


2 

n 


fT 

0 


rT 

h(x) h(z) R^(z-x) dzdx 

0 


2.3-4 


The optimum signal-to-noise ratio can be found by minimizing the 
Lagrangian: 


L = 


T T 

I I h(x) h(z) R^(z-x)dzdx - p 
0 0 



2.3-5 


The resulting variation yields: 

fT 

ho ^ s(T-x) 

0 


2.3-6 


The filter which satisfies this relation will maximize the signal-to- 
noise ratio for a known signal in any noise with autocovariance function 
R^(x). Equation 2.3-6 is, of course, a Fredholm integral equation of the 
first kind which is solvable only for a restricted group of covariance 
functions hov^ever, we can adequately approximate the 

integration by replacing the 0 to T limits with -» to +“, then the 
Fourier transform of 2.3-6 gives immediately 


- S„(s) 


where s = iw and S^{s) is the actual power spectral density function 
of the noise. This matched filter is then just the white-noise 
matched filter co nvolved with the actual noise spectrum. This result 
was derived for the limit T->±», but a detailed derivation shows that it 
holds wherever the data samples are uncorrelated, which was determined 
from the zeros of the noise autocovariance function. 
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As In the previous section, a best estimate of the signal amplitude 
exists in the presence of non-white noise. That estimate is given 
by: 


A = 


fT 

, Qh(x)r(T)dT 

h(T)s(T)d,T 

,0 , 


where h(x) is the solution of: 


s 



h{x) dx 


2;3-7 


2.3-8 


Comparing this result with 2.3-6, we see that the optimal whitening 
filter is the best v/eighting function for the correlation detector and 
the amplitude estimate in the presence of non-white noise. 


2.4 Karhunen-Loeve. Filtering 

The emphasis in the preceding section was an additive white noise. 

Since this is often invalid, v/e derived a test based on the noise auto- 
covariance function. For the v/hite noise case we considered a flat 
band! imi ted spectrum and found that appropriate uniformly spaced amplitude 
samples were statistically independent. For colored noise we considered 
the continuous sampling limit and wrote the detection equations as 
integral relationships. However, uniformly spaced samples in colored 
noise are correlated and the sampled case is difficult to evaluate 
explicitly. There is, however, another method which can be used to 
generate statistical ly independent samples. While these are not 
amplitude samples, they can be used to construct the same detection 
and performance tests as previously described. The approach ^used will 
be to expand the signal in a series of functions which are orthogonal 
over the region, 0 to T. 

The functions v/e seek are a set of f^-(t) 's with the normality 
.condition: 
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which is identical to 2.2-4 except that has been replaced by 
the eigenvalues and the signal samples have been transformed by a 
weighting function similar to the v/hitening filter of section 2.3. In 
this case, however, equations 2.4-2 through 2.4-4 can be written as 
sums over the time sampled values with no loss of generality, hence, 
with no degradation In performance caused by correlated samples. 
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2.5 Digital Filtering 

The -transformations of sections 2.2 through 2.4 can be written as filter 
transfer functions. ■ Additionally, empirical methods can be used to 
synthesize a desired transfer function and the equations of those 
sections can then be used to evaluate the error rate and detection 
performance. This latter course is often followed ^hen the sampling 
rate is constrained by some considerations other than those requiring 
uncorrelated noise samples. Typically thedesired filtering is 
matched to the sample rate and the signal dwell time by the Nyquist 
theorem and we wish to evaluate the detection performance of such 
systems. Additionally, it may be desirable to further filter the 
data to improve the signal-to-noise ratio based on the observed noise 
spectrum. In this section we will discuss how such a transfer function 
could be synthesized and then derive the algorithm for converting that 
analog transfer function to a digital difference equation. 

Given an analog impulse response function H(S), the difference equation 
for the filter function can be derived. Also, given the nominal 
characteristics desired, the transfer function can be synthesized. 

Both of these techniques are described below. 

The frequency response can generally be described as a series of 
first-order filters. The transfer function of a low-pass filter is: 

= 2.5-1 

v/here a.. = 2irf.. , is the gain of the filter. 

f.j = the corner frequency of the filter (Hz). 




That is, the response of a low-pass filter is flat for f < f ^ , and falls 
at 6 dB per octave (linearly on a log A-log f graph). 

For a high-pass filter, 

"" ^iTTT^ 2.5-2 



and the slope is the same as before. Higher order forms of these 
filters have transfer functions which are powers of the above H(S)'s, 
with the exponent n equal to the order of the filter. That is, a 
3rd order high-pass filter is: 

H.(S) = 

and its response slope increases by a factor of n (3 in this example): 



Finally, a circuit which can be described by a series of such, filters 
has a transfer function which is a product of the elemental H. (S) 
terms, and a response curve which is a series of line segments with 
n(+ 6 dB) quantum slope changes at each characteristic frequency. 
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The definition of the corner frequency is important here. If we have 

i 

defined the f^ points as the traditional 3’dB or half-power response 
frequencies, then the values used in 2.5-1 and 2.5-2 for a^. must 
be altered *somev;hat. Note,, for example, that for a single order 
filter (m=l), v;e have the half-power response at 


H(S) h'IS) = i ^ S=a. 

as expected. For an order filter, we find 


2.5-5 


HCS) hIs) = 


1 

2 


= { 


1 

at + 
1 




^ $’■= - 1) a^ 


2.5-6 


and similarly, for a high-pass filter: 



2.5-7 


Given the transfer function H(S), the difference equation can now be 
determined as follows. First, transform the frequencies. Since we 
desire the digital equivalent frequency, determine by: 

2irf.T a-T 

A. = TAN i—2^) = TAN (-^) 2.5-8 

where T is the sampling interval, = 1/SR {SR is the sample rate in 
Hz). Second, transform the H(S) function to an H(Z) function by 
the substitution: 

2.5-9 


This transformation preserves the square-magnitude response of the 
system except for a warping of the frequency scale as given by the 
first relationship (2.5-5). The advantage of this transformation is 
that aliasing is not introduced by the digital representation, thus 
avoiding the necessity of "guard" filters (anti-aliasing) which would 
result in a digital filter of higher order than the original analog 
transfer function. 
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The technique thefi is to express H(Z) in terms of a polynomial in Z, 
as: 


H(Z) = G * ^ 


with the order of Q{Z) equal to or greater than P{Z). Most easily, 
the substitution used is: 


,A, + 1)(Z + U.) 

(A: + s) - (z;;, 

(A. - 1) 

where and = A. + 1 


2.5-11 


The numerator and denominator of H(Z) are then divided by Z*^, v/here 
n is the order of the denominator, resulting in a transfer function 
which is a ratio of two polynomials of equal order in powers of Z~^ . 

Equation 2.5-9 is used for S terms in the numerator. The general form 
of this H(-Z) is: 

Z-n(2+i)n-m ^ 

H(Z}= K G. 2.5-12 

(1 + I s. z"’) n p. 

i=l ^ i=l ^ 

where : 

. A^- is as defined previously; 

n * the number of elemental filters; 

m = the number of high-pass filters, and the A. 's are ordered 

.V 

with the low-pass frequencies first. 


and 


are the expanded coefficients of the product function: 
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n 

n 

i=l 


(U. + Z) 


such that: 

n 

Si = I U4 

’ ,1=1 ’ 


n 0=1 

Sg = I u. ( I U.) 

^ . 0=2 ^ i=l ^ 

n ■ k-1 o*-l 
S3 = I Uj^ n u. (I u.)] 

0=3 0=2 ^ i=l ^ 


n m-l £=1 


j-1 


S^ = I U I I [.... ( I A )] 

^ ra=r £=r-l k=4-2 i=l 


2.5-13 


Again from our example 2.5-4, the H(S) transforms to: 

„(Z) = . p (Z- . 1)( , ?±1 . L^ S 

Z'^(Z+U^)2(Z+U2)(Z+U3) 


•‘2.5-14 


The Z-transform corresponds exactly v/ith 2.5-12 if we consider the 

four elements of the filter as having frequencies (in' the. original 

form) of 1, 4; 40, and 40 Hz. That is, U.j is repeated, but treated 

as if it were two different terms. Ignoring the constant factor 
/ 

temporarily: 


H(Z) = 


1-2Z"^ - 2Z"^ - Z"^ 
4 1 

1 + I Z'l 

i=l ^ 


2.5-15 
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where S-| = 2U-j + U 2 + 


$2 = + 2U^U2 + 2U^U3 + U2U3 

S 3 = U^U^U 2 + U^U^U 3 2 U^U 2 U 3 


2.5-16 


Now that we have the Z-transform of the H(S) transfer function, the 

-1 

difference equation for the system can be written. Noting that Z 
is the unit delay function, and writing H(Z) as: 


1 + l.j. 

H(Z) = - 

1 + I s. Z"' 

i=l ^ 


- Y(Z) 

■ m 


2.5-17 


where Y(Z) is the Z-transform of the output, and X(Z) is the input 
transform. Inverting the transform we find: 




n I i 

. + y T. X. . - y s. Y. . 


2.5-18 


where Y. is the jth sample of the output series, and X- is the 

J J 

corresponding jth input value. Continuing our example. 


VXn-2\-r2X„-3-Xn-4-Si''„-rS2Xn-2-S3Xn-3-Vn-4 


2.5-19 


It is interesting to note that because the transformation 2.5-9 is 
bilinear, the difference equation v/ill always be the same order in 
X^_-j 3nd terms (except for canceling of some Z^ terms by the 
expansion of (Z+1 )”’(Z-1 )^). 


To find the T.. coefficients of equations 2.5-17 and 2.5-18, we must 
expand 

(1+Z)"""’(1-Z)"’ = 1+T^ Z+... + T^Z" 
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2.5-20 



but the gain cohstant in 2.2-12 of 


m n ^ 6. 
n A. n (“) = COEFF 
i=l ^ i=l M . 


2.5-21 


must be included as a factor of all the'X^ terms. 

Expanding 2.5-20 can.be done using the binomial expansion subroutine 
attached to expand (l+Z)^"*^ and (1-Z)"^, and the polynomial product 
routine to find the resulting terms. Then if we redefine the sub- 
script of T by one to absorb the 1 on the r.h.s. of 2.5-20, we can 
include the factor 2.5-21 easily into the definition of : 

(1+Z)‘^“"’(1-Z)"’ - + T 2 Z+... + Z" 2.5-22 

if we also redefine the and set S-j = 0, then the relations 
2.2-17 and 2.2-18 can be written: 





2.5-23 


^0-i+l 


-This redefinition of S.^ has been included in the attached algorithm 
(see section 3.1). 
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Section 3 of the first report covered in detail the basic survey 
point-source processing scheme. However, in order to successfully 
complete the sky survey and to define the detection test gates, a 
number of peripheral routines are needed. This section will discuss 
several of the most important of these routines. Some routines are 
simple, such as coordinate transformations used for positional 
matching. Others, although complex and worthy of discussion, are 
very specifically written for each mission. These routines generally 
are part of the primary data processing system and are assembled from 
the formula and algorithms of section 2 of this report and section 3 of 
the previous paper. Yet, other routines which are a. part of the data 
base merging are decisions and tests for specific types of astronomical 
sources and depend on the sensor bandpasses and sensitivities 
and on the spectral characteristics of the sources being searched for. 
Some routines which are only peripheral to the primary detection 
scheme are so basic and important that they are worthy of individual 
discussion here. The complete set of programs designed and tested 
on the Celestial Mapping Program (CMP) data will be published at 
a later date when that task is completed. 

Two basic programs will be covered here and one front-end detection 
scheme used on a previous survey program. First, we will discuss an 
application of the digital filter design scheme of section 2.5. Then 
the calibration of infrared detectors is discussed and a routine to 
evaluate the spectral response of a detector plus filter combination' 
to a variety of stellar spectra. 

3.1 A Digital Filter Design Aid 

The discussion of section 2.5 covered the algorithm for digital . 
filter design; here, we consider the specific use of the following 
routines. The program attached does two things. -.First, the 
coefficients of the difference equation 
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N+1 
J i=l 






where the X. 's are the input data, Y^'s the output data, and 

and the filter coefficients determined from the desired frequency 

response. 

Second, the routine creates a pair of sample response sequences. 

One is the impulse respdnse function of the filter. If the filter 
characteristics were chosen to duplicate the response of a 
detector and its electronics, then this impulse response will model 
the radiation-hit' response. The other response is the system 
reaction to a square wave. Since the duration of the square v/ave 
is equal to the point-source dv;ell time, the response is approximately 
the same as a source signature. 
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PROGRAM TRNSFN 


76/76 OPTsl 


FTN 4,5 + 4 H 


05 / 06/77 


X program trnsfnunput, output J 

COMPLEX C^SUM,COEFp,T 
C 


5 C , A*AA 

C • THE FOLLOWING DIMENSION STATEMENT IS USED TO CREATE — 

C A sample, response sequence for a rectangular input aaaa 

c equal in length to the star dwell time, and the impulse **** 

c response.' sequence, ■ 

I 0 ' C AAA* 

r dimension XX( 200 ),XPl 200 ),YY(aoO),YPC 200 ),TIM( 200 ) ''*** 

C ■ ' AAAA 


Caaaaaa*aaaaaaaaaaaaaaaaaaaa*a*aaa*aaaaa 

■ c • ■ 

15 C the routine is dimensioned for transfer functions of total order 

C 9 OR LESS, FOR HIGHER ORDERS, CHANGE THE D.IMENSIONS .OF ALL THE 

C following TO NtJ, WHERE N IS THE DESIRED ORDER, • 

C 

dimension CnO),SUM( 10 ),TClO)f xaO),yUO),<i( 10 JfRS( 10 ),TRnO) 

20 C 

PRINT 99 

99 FORMATCiHRJ 

AlslHC a A 2 = 1 HS $ A 3 = 1 HT $ A 4 =IHF 
C 

25 C N AND M. DESCRIBE THE ORDER oF THE TRANSFER FUNCTfON, 

C N s the total order 

C M = THE OROEH OF THE H 16 H«PASS FUNCTIONS 

C 

c the program will loop for new transfer Functions, to end, 

30 C set M » 0 , (THE LAST DATA CARO CAN BE A BLANK TO STOP) 

C 

• 50 CONTINUE 

• READ 100 , N,M 

100 pqRMAT ( 215 ) . 

35 IF(N',EQ, 0 ) GO TO 51 

C t 

C the sample rate SR governs THE FREQUENCY WARPING FOR THE S TO Z 
c ■ transformation, 

C THE DWELL TIME IS DETERMINED FROM THE SCaN RATE SCNRTE AND 

40 C THE DETECTOR SIZE SIZE, 

C , • 

C SIZE And SCNRTE ARE VARIABLES USED ONLY FOR GENERATING 

c the test cases, 

, C 

A 5 read iQf, 5 R,STZE/.SCNRTL 


M 

C^ 



n o o o o o o 


■ l,N?lNTC(S«ttslZE/SCNRTE)+Q,S) 

H5T the input parameters 
PRINT 20WN,M; SR, UN 

201 FORMATUHI ,//,20X,*N s M = SAMPLE KATE - *,F 

nOe«/* DWELL = 

the corner frequencies for the filter elements can be specified 
AS complex (REAL + IMAGINARY) 

THEY should be ORDERED WITH THE M HIGH-PASS ELEMENTS LAST 


ts> 



PROGRAM TRNSFN 


76/76 OPT=l 


PTN 


05/06/7"? J7t01,5’i 


60 


65 


70 


75 


80 


85 


90 


95 


100 


C EACH ELEMENT HAS A (COMPLE^OGAIN WHICH CANNOT 8E ZERO 
C THE gain can be SET AS (1,00,0,00) 

c the orper must be egual to or greater than one 

C • , ■ 

C THE FREQUENCY CAROS ARE FORMATTED (^FtO,6,J5) AS FOLLOWS; 

C FREO(REAL)rFREQCIMAG),GAIN(REAL),6AlN(IMA,G3, ORDER 

C 

PIS3, Iill59a653579 

COEFF-CMPLXn, 0,0,0) 

■ NT»0' 

00 i 1=1, N 

READ 102, WDR,WpI,GAlNR#GAINI;NORD 
102 F0RMAT(4|F10,6,I5) 

PRINT 205, I,WDR,W0I,NDR0 

•205 FORHATHH , 10X,*WD(A, 12,*) = * , 2F 1 0 ,6, SX , *0RDER .= y^,J5) 

NN=NOROv 

IF( (NTffM) ,GE.N) go TO 7 
KDR=WDK/(2,**C 1 ,/NORD)«l .) 

W0I = W0l/(2,**(i b/N0R.D)^1,) 

GO TO 8 

. 7 W0R = ‘rtDRA(2,*AU,/N0RD)<-l,) 

KDI = WDlA(2,A*(i,/NORD)-.i,) 

8 CONTINUE 

KA.I = TAN(hDl*PI/ SR) 

WAR = TAN(^^DR*PI/ SR) 

C - 

c •, list the frequencies 

C 

•131 CONTINUE 
NT=NT+i 

PRINT 202, NT,WDR,WDI,I,WAR,WAI,NT,GAINR,GAINI 
202 FORNATUH , 1 0 X , *WO ( * , 12, * ) = * , 2F 1 0 ,6, 5X , *WA ( *, 1 2, iv ) = *,2F10.,6,5X 
1 , aGAINC*, 12,*) ? A,2FiO,6) 
C(NT)=CmPLX(WAR».1,,WAI)/CmPLX(WAR'+ 1,,WAI) 
C0EFF=C0EFF*CMPLX(GAINR,GAINI)/CMPLX(WAR+1 ,,WAi) 

IF(NT,LE, CN-M)) .COEFF=CO£FF*CHPLX(WAR,WAn 
NN = NN**i 

IF(NN,GT,0) GO TO 131 

if(nt,eq,n) go to l«) 

• 1 CONTINUE 

HI CONTINUE 

C . 

C . LIST THE C(I) terms and THE VALUE OF COEFF 
C 

TRrC 


to 

00 



o o 


uo 


r 




ao'i FORMAKiHO,//, aiOX#Al,*(*y.I^^l B ) 

OeTFRMINfc- THE S(J) COEfFlCIEMTS OF THE Y(J-U TERMS 


DO 2 iBliN 

2’ SUMa)sCMPLX(0»0,0,0) 

DO fl J=1#N 
suFcUBSuMcn+ccn 
lFa,EQ,N) 50 TO U 
DO 3 0*2, N 

IFCIW-^I. ,Lfc,N? SUM(0)BSUM(J)+8UMN'-U^CatJ^n 


tv) 

to 


ORIGINAL PAGE fj 
OF POOR QUALHa 


Trogram trnsfn 


76/76 OPTsl 


fTN 4s5-fr4ta. 


OS/06/77 U.OUSi 


ils' 


UO 


1^5 


130 


135 ' 

li'ao 


1«5 

150 


155 


.3 CONTINUE 
a CONTINUE 
00 5 I»1,N 

■ 5 SUMCN + 2-U?SUM(Nf 1«I) 

SUf<(n:5CMPUX(O,O,Oj0) 

c 

c ■ ,uisT the sen coefficients 

C , ■ 

NN?Ntl 

PRINT 203, (A2,nSUMU),I = l,NN) 

C 

C ■ EXPANO the NUMERATOR TO DETERMINE THE TU? COEFFICIENTS UF X(J«U 
C 

C note THAT TCI) HAS EXACTUY (N’M+ 1) + « N+1 TERMS 
C ' • ' 

NPfN'wM 

. CAU BlNEXPCXfNP, + nOO) 
lOIMXsNP 
NP5M 

CAUl BINEXP(y,NP,>-nOO) 

IOIMY=NP 

CAU PMPY(X,IOiMx,YnOlMY,Z,IOIMZ) 

DO 6 nniDIMZ I 

T( naZ(I)*CO£FF 
TRCn=REAL(T(I)) 

RSensREAUSUMdJ) 

6 continue 
C 

C NOW LIST THE T(U COEFFICIENTS OF THE YCJ-l) TtRHS 

c • 

PRINT 203, (A3,I,T(nrI-niDlM/:) 

C 

C NOW AUU THE FACTORS ARE DETERMINED, 

C Wf may now calculate THE RESPONSE TO AN ARBITRARY INPUT SEQUENCE 
C OF X(I) VALUES AS ^ 

C ' ' ■ 

C YtJ) =» SUM(I = l TO NtU OF (( TCn’*X<J-»I + l)"SCn*YCvJ»'I + l) )) 

C 

C , ■ 

CA**A*A*A^*AAA****A****A**A*AA:^.AA*kAA A*** 
C • A*** 

C the following CAROS -HARKED aaaa CREATE A SAMPLE RESPONSE aa** 

C -■ SEQUENCE, THE INPOT RECTANGULAR PULSE IS CREAJED IN XXCI) *a*a 

C' WITH A LENGTH EQUAL TO THE STAR DWELL TIME FOR THE INSTRUMENT Aa** 

f' T^PIIT TMF PV SPOmKS- SPf.llPMrr T<; TM VVfT*l Tut- TMPHISt- ***'*. 


OUTPUT IS NUmEJERLD SEQUEh '{iLLY AND THE SAMPUE TIME IS **** 

GIVEN IN MILU8EC0ND6, 

DO u isi,aoo 

TIMd)sO,0 

iFd.GEdO) TIMd) = lOOO«*d«ao)/ SR 

xx(n=o,oo $ xp(i)=o,oo **** 

YYd)=;OsOO $ YPdJaOjOO 

CONTINUE ***A 

LL=19+LN 

DO 12 Is20fLL 


■PROQRAM TRNSFN 


• 76/76 ' OPT51 


FTN ^ 95 +«i 4 


, 05/06/ff J.7,01,5I$ 


i75 


lao 


J85 


no 


195 


200 


205 


2 i 0 ‘ 


12 xxaj*i,oo 

. C ■ 

c - ■ ■ 

C generate GAUSSIAN NPISE 'AND ADD IT TO THE INPUT DArAj 

C • ' 

CALL RANSET(6632U? 

A T 0 9 

NOISES I S SNRslOO* 

NOISEnl $ SNRsiO, 

NOISE=0 

IF(N0I$E',EQ,0) GO TO 123 
DO 122 isio^aoo 
GN0ISE=0*0 
DO 121 J = W12 

121 GNOISE = GNOISE+ RANpCA) 

GNOISE k (GNOISE!^6sO)/8NR 

xxcnsxxnitGNoisE 

122 CONTINUE 

123 CONTINUE 


XP(lO)-lpOO 

iLsNU 

jFsN-f 2 

YYKst9999, SiYPMs»9999, 


C 

C 

C 


U 


14 


15 


300 

111 

A A 


DO H JaJF^aOO 
DO 13 IsWiL 

YYCJJsYYC J) + TRCI) axX(J-I + 1)-.RSU)*YY(J-»I*1) 
YP(J1?YPU)fTR(I J*XP(vJ-^I + 1)-RSCI)aYP(J«U11 
CONTINUE ■ 

YPNsYP(J) 

YYMsYY(J) 


IF(YP(J) jGT,YPM) 

IF(YY(J) ,GT «YYNO 
CONTINUE 
DO 15 Jsi,200 
YY 1 JJsYY(J)/YYN 
CONTINUE 

PRINT 300» YYM^YPM , 

FORPaTUHI , 10X,*YYMAX = a,F10,6,a y'pmAX s a , F 1 0 , 6, // ) ' ' 

PRINT lllf (I»TINCUfXX(n,YY(I),XP(I),YP(I)^I = l,200) 


$ YPCJ)=yP(J)/YPM 


1 


FORMaTUH 

xp 


I NSEC X 

YPa,/, (( 3X,I3,5(5X^FiO,3)))) 




AA* A 
A A AA 


A AAA 
A A A A 
AAA A 

AA AA 
A AAA 
A A A A 
A A A A 
A A A A 
A A A A 
A A A A 
A A A A 
AA A A 
A A A A 
' A A A A 
A A * A 
A AAA 
A A A A 
A A A A 
A A A A 
A AAA 
A V A A A A A 


Ol 

ro . 


c ) 
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GO TO 50- 


SYmBOUC reference’ f>'AP (Ran. 

ENTRY points 

^5 trnsfn- 



3.2 Infrared Filter Calibration 


The calibration of infrared brightnesses is the single most 
difficult aspect of a scanning sky survey. For ground-based 
point-and-integrate systems, it is possible, and in fact common, 
to make all measurements of the same signal-to-noise ratio by varying 
the integration. Since the amplitude uncertainty was shov/n to be 
a function of the signal-to-noise value in section 2.2, it is clear 
that uniform photometric accuracy is readily achieved. For 
survey instruments, the uncertainty of the initial measurements is 
inversely related to the signal-to-noise value, applying a fundamental 
limit to the accuracy of the survey measurements which varies both 
with brightness of the source and its location in the sky. 

Further complicating the problem is the fact that the sources have 

a wide variety of spectrum so that the broad band detectors typical 

of infrared survey instruments do not have a well-defined intrinsic 

calibration. It is possible to calibrate the detector voltage 

2 

in terms of the watts per cm it receives. However, if the survey 
is measure sources in several colors, or if the calibrations are to 
be derived from measurements made in a difficult wavelength region, 
the measurements must be referred to a spectral intensity. The 
wavelength bandwidth that is needed, however, is dependent on the 
spectrum of the source being measured. Furthermore, the effective 
wavelength of that measurement varies with the input spectrum. 

The units of the brightness measurement are another problem. The 
most useful form is the brightness magnitude, defined by 

m = -2.5 log,Q (/) ~ 3.2 

where B is the observed brightness and is the reference value. 

This reference is different for every filter, since it is defined as 
the response that filter-detector would observe from a particular 
"standard" star - the archetype is a-Lyra, which is defined as a 
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10,000°K black'body source with an angular diameter of 1.5697 E-16 
steradians. The great benefit of this magnitude measurement is 
that vie skirt the question of effective bandwidth. These magnitude 
measurements still need an effective wavelength, but for blackbody 
spectra at wavelengths less than 50 micrometers, the effective 
wavelength varies only very slowly until the source temperature 
falls belov/ 500°K. Finally, the magnitude measurements defined 
by 3.2-1 can be used inversely to find the equivalent blackbody 
color temperature if measurements are available in tv/o or more 
hands. 

Figures 3, 4, and 5 show the variations in bandwidth, effective 
wavelength, and magnitude difference for three infrared filters 
similar to ones commonly used in previous surveys and measurements. 
The results were derived from the attached filter calibration 
routine which is self-explanatory. 
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iACCOUNTf SPTHEM, T302S 
I HAP. CP K, 
iFTN'. 

iREKINO. PUNCH, 
iCOr-Y, PUNCH, OUTPUT, 
■EXIT, 


C 

c. 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PH06RAP IN TORT CIKPIJT , OUTPUT, PUNCH , TAPEbs I NPUJ . TAPE6 = 0U T PUT ) 

PLANCK2 (HCDIFIEC) 

IKTEGPATIOU can be done WITH OR wl rUDUT THE PLaNCK rUNCTlOH* 

PHCURaH PUR INTEGRATING THE PLANCK ULACKHODV RADIATION FUNCTION 
OVER an Interval CETERMIHEP and attenuated by FlLTEH-SENsniVlTy 
TYPE FUNCTIONS, 

data sHoUl.0 BE IN order OF INCREASING WAVELENGTH. 


c 

NFUNCT 


NO. 

c 

IRITF. 

Z 

1 - 

c 



0 S' 

c 

NF 1 UX 

z 

1 

c 


z 

0 

c 

IRAVTR 

z 

1 « 

c , 


m 

6 “ 

c 

NUphAL 

z 

1 * 

c 


z 

0 *“ 

c 

IPLGT 


1 ft 

c 


z 

0 CO 

c 

ntrans 

z 

NC. 

c 

H 

z 

NO, 

c- 

N 

z 

NO, 

c 

L 

z 

NO, 

c 

EXPEP 

z 

N A F t 

c 

TYPE 


IPEN 


DO HOT WHITE, 


intensity, radiance, PLANCK INTENSITY, 


PL/rCK FUNCTION NOT TO BE CALCULATED, 

TRAH3 PISSIUN data DECK CONSISTS OF WAV, AND TRANS, 
TRAf«'GMlSS ION DATA DECK CONSISTS OF TRANS, DATA ONLY, 
KCPfALnEU FUt'fCTJOMd) * FUNC.lION^a) 

CO f'CT NCPhaLUE 


OF SETS CF transmission DATA PER RESPONSE FUNCTION, 


Card i 
CAPO 2 
Card 3 
card ii 
CARO S 
CARC 6 


CF TFMF. TO BE CaLCULA(F.I) (IF rMERE ARE NONE Ls 1) 
OF EXPFRlHEMr DATA, ETC, FURhaT ^A6, 

TIFYING KAKt C-F DATA. FORMAT 4IA6, 

N' F U N C T , I R 1 T E , N F L U X , I W A V T R , N 0 H M A L , I P 1, 0 r 
N'T RAMS 
P f , I , 

U^'P (isi.L) 

E>PfR (|<sl,<U " hAV,RE 8 (Jsl,,N) 


(6131 
(Hi) 

an ,213) 

(9FS,2) 
(i|A6/a2F6,9) ) 


T'fPf = « ('.Av,TPAN oh TRaN (1 = 1, N) ( AS/ t J pf 6 , ) 


* * :A ^ 


■A ii -k -A 


A 'A AAA A, A A A AA A AAAAA A A AAA 


reference flux Is the INTEGRAIEU FLUC of a 1 0,000 DEGhtt. 6,B, OF 
Sl/E l,b6S7E-16 SIR (ALPUa LYRA) . , . THF /CRO MA(jNX fUOE REF 
and is integral of Ei (LAh'BuA f 1 0, OOO )*R (LAMBDA ) *>PLAMUUA*0HEGA 




hj Q 

§1 




6q 


O'? 



o o o n o o noor^noooooooooooooooooo 


OEFJNtO wavelength IS THfc FLAT KESF’ONSE , EFFECT I VE WAVELENGTH 
BDWDTH 2ER0 IS THE FLAT RESPONSE SANDWIuTH. 

Wavelength is the true effective wavelength for b,b» spectra at 
' the given tehferature 

the true effective wavelength is the integral of H(LAMDPAfT)^ 
R(LaMHCA)*LaH20AACLaMBDA divided fiY THE INTEGRAL OP D(LAHt3DAfT) 
*K(LAf'BDAI*DLA.vDDA 

BDWDTH IS THE TRUE BANDWIDTH- AT THE TRUE EFFECTIVE WAVELENGTH FOR 

the given tehferature 

The true effective bandwidth is the integral of B(LAHBOAfT)* 

■ R(LAF'bDA) aDLA:-'PDA DIVIDED BY BCLAH-BDA EFF,,T) 

/ BDwDThDF IS THE PROPER BANDWIDTH FOR THE DEFINFD EFFECTIVE WAVELENGTH 

BR HAG IS -*2,b*LOC(0(LAHODA EFF / T ) /BO (LANuDa 0,10)) WHERE BCLAM0OA,T 
IS THE B,E, EF I SSI ON- 

. COL' HAG IS ~2,5ALCC( JNT FLUX/REF ERE WCE FLUX) 

•. • INT Fl.U}( IS THE INtcGRATEO F.lUX ON THE DETECTOR IN W *CH^2*STR« J • 

VIEW flux is the integrated flux multiplied by THE FIELD OF VIEW 

= ii,bE-06 SIR AND IS IN w*CM^2 

( 

REAL MAGFaC 

DIMENSION WAV( 150) , RES (ISO) , TRAiJ(lbO) , TEMP (ISO) , W.( ISO) ,wX( TiO) , WXX 
HISO) ,«ESP (ISO) ,WWX( 150) f TYPE (ii ) , E XPEH (a) 

NJ = 0 

OHEGA=l ,Sb97E-16 
, PI = :S1^I592E5^ 

ICOL=0 

L ^ 0 i 

READ (5,9i|) NFUNCT, IRITE,NFLUX, lAAVTR, NORMAL, IPLQT 
9^ F0RMAT(6Ii) 

95 NvJ = Nj4l 

NFLUXsO 
DO 92 1 = 1, /i 
92 TyPEU)=SH ' 

LOLD=L 

IC0L=IC0L+1 



Nf^sO 

^tAD (5,99) NIRAN'S 

99 FUHMAr(ll3) - 
HEAD (5,lUO) 

100 FuRi"Ar(Uiraii) 

( 1 ) = 1 00 0 0 , 

c Insert te<<p oEKifv'iRG cards? e=xx, temp etc, 

l = 70 

TE>-(P(2)s100. 

DO ?0 1=5,412 

20 Ttt^P(I)sTl.'’'P(I'*l)4lO» 

OC 21 1=93,52 

21 lLHH(l)=TtMP(I-i) tSOt 
00 22 I=53i70 

22 ■ TtPP (I )=7£PPn-n+5C0, 

Hi=l ® N2=R 

READ (5^ 10 1) (EXPER(K)ft<5l»9)^ (WAV U ) » RES (I ) , JsNi , N2) 
101 f]0KMATC9A6/(l2F6,4n) 

l!l=R + l 

DO' 500 I»L1, 150 
HAV(I)=0.00 
RES(I)::0,00 
1HANCI)=0,00 
RfcSP(n = 0,00 
500 CONUNUE 
. N2 = N 

9e NN=NN+1 

IF(i'i-l) l^fi 

3 If UWaVTR.EG, I ? GO TO 601 

read (5,112) (TypE(«),K = l/9Jf (TRAN(I),I=:l,N) 

112 FURHAT (4lA6/(l2F6,4j) ) . 

GO TO 600 

601 READ (5,112) (Typ£(K),K=W9)/ CWAV(I), rRANU),I = lfN) 
600 CONTINUE 

DO 2 I = 1 

2 HESP(I)=RES(I)*TRAN(I) ' 

GO TO 11 •» 

I DC- 9 1:^1 I 

TKAiv(i)=0, 

9 HESP(I)=RES(I) 

687 CUlvTIf'UE 
U Cun'IInUE 

lF(NDRf^Al,FG. 1 ) GO TO 700 
urGHR=l. 

M = 1 


4^ 


o o 






h'SPACt = G 
L i N fc. •» 0 ■ 

WHUF (6f212) CEXPFR(K),I^ = J CTYPE(K),Kj:1/^) 

•'212 FORf- A r ()H} ,37Xf (|A6/ '4 1 X f 1 9HT ABUE OF IR'PUf DATA/ 3flX, 0A6//20X, 

UOHi- AvEL£NC-rH,3XK 1 lHru.VCTlO^(U/3X. UH>'UWCT10N(2) fSX / jONWAVtlENGTH 
2, iX, 1 1HFU!\CUC^J( n ,3X, UHrURC I IUN(2)//3 
I&G = N’l i 

> DO lU X»lRG,r4D 
NSP/4CE=IJ5PACE+1 
LIK;fc=:LlK£-H 

IF (b"t-SpACr)20'T,206, 206 

206 WRITE (6^210) Wa V ( I) , RES ( I) / TR Afl ( I ) # WA V (I + 3^) ^ RES Cl + 3^3 , TRAN (I + iS) 

210 FURhat(22x,FS,2,7X,F7,^,7x,P 7,^, lOX»Fb,2,?XfF7,^,7XrF7,«) 

GO TO 

20^ NSPACF.= 1. 

IF(i5-LlNb)2l9,2H,2H 

ZIU WHITE (6/211) WAVU?;ReS(I)/TKAiMa)/WAV(I + 3b),RtS(If3‘i.)/TRANCI + 35) 

211 FORi^AT(lHO,2lX/r'3,2,7X,p7,^,7X,F7,^,lOX/Fb,2,7X,F/,A|,7X,F7f'il) 

GO TO 

219 Nl=l+LIN£*1 

If (N2-M )250,220,220 
n-T CONTINUE 
GD TO 250 
700 BiGHPsRESP(l) 

DU 40 1=2,K 

IF (BIGHF-hESPa) )32, 30,30 

32 • 8IGRP = RESP(.U 
30 CONTINUE 

00 33 1=1, N 
RESP(1)=R£SP(I)/0IGR.P 

33 CUNUNUE 

■ ' Nl=i 

720 NSPaCE=0 
■ UN£ = 0 

•WKiTE (6,712) (EXPER(K),K=i,iU, (TYPE(K)fK = l,iO 
• 712 FORMA f (IHI ,/16X,AA6/ SOX,19HTAOL£ OF INPUT DATA/ «7X, 0A6//10X, 

1 10HnAvELENGTH, 3X, nPEUNCTlON'(l),3X, 1 IHFUNCT10 n(2),3X, 

21 IHFUmCT , ( 1X2) ,SX, ICHWAVELFRGTH, 3X, 1 Ihf UNCTION (1) r3X, 

. 31 U-F UNCTION C2) ,3X, 1 IMHUvCT , (1/2) /SIX, 1 OhNOKHAL I /.EO , 47x, 

^IOHNGPMaLIZEO//) 
s lBG=Nl $ IM)=N2 
/ 00 6iq 1=IPG/INC 
NSPAr.F=NSPACE-U 
• UNF=LlKE-fl 

IF (b-NSPACF)70a, 706,7.06 



706 i-^HUE, (6,7103 V ( I ) , RES (1 3 , TKAn ( 13 , R£SP ( I ) ? WAV ( I + 35 ) , RES U +3b) , , 
iTKAN-n + 35 3 ,KE5P(I + 35) 

710 FUR«AT(12X,F5.2,3C7X,F7,^-^), i0X,F5,2,3(7X,F7,^)3 
GO' TO' 6li( 

70^1 N’SPACt^i 

lF(35~UW£37l9,7i^,7H 

71^ WRITE (6,71 n 'l^'AV(l3,RESn),TRANa3,RfcSPU),WAV(I + 35),RfS(I + 3S), 

1 TRAlV(I + 35) ,RESP(If3S) 

711 PORhAT(iH0,nx,F'5.?,3(7X,F7,«),10X,P5,2,3(7X,P7,A)3 
GO TO 61/i| 

719 Kl = I + UNE-l 

U(NP-M) 750,720,720 
61/4 COnTIHJE 

7S0 continue 

nil PORr*AT(F10,7,Fl0,83 
250 IF(nFLUX,CG,03 GO TG 307 

wKJ rE(6, 1 07 3 nEPP(I) ,l5l,U3 
IQ7 FOKt'iATC// 55X, ilHTC»«‘.Pt‘RATURe//(9Fl3,2)3 
600 J=0 

6 ‘ J=J +1 

IF(NrLUX.EG,0) GO TO 307 
I3U ^ 1 = 1,.N',1 

4 CALL 5i^CTK« C TEf-F? ( J 3 , WA V ( 1 ) , RESP (1 3 , W U 3 , WX ( 1 3 , WX X ( 1 3 , Z INT 3 

GO TO iOl ; ■ 

307 GO- 302 1 = 1, N I 

wXC.i}'='RESFa3 ■ 

A>rX(l-3f:wX(l3 /PI 
302 OOn'TTnUe 

• • IFCIRITE.EO.O) Go To 300 

':."301 NSPACC = 5 
UINE-95 
00 2)6 1=1 ,N 
NSPACE=nSPACE+1 
UIsE'3LINP + 1 

IF C5-HSPACF 32^10,2^2, 2/J2 
2A2 wRITF ( 6 , 1203 W AV ( 1 3 r K’X ( I ) , NXX ( I) / W ( I ) ’ 

120 FORN4T (30X, l PE 1 2 . 5 , /J 1 PF i 2 , 5 , 5X , 1 PE 1 2 , b , 9X , 1 PE 12 , 5 } 

GO TO 2V6 
2/lG N'SPACt = l 

lK^b-l.lGE32«/4,2/46/2<46 
2/|6 isRlTE (6, 1 223 N'A V ( 1 3 , XX ( 1 3 , V.XX U ) , 'W ( J3 

122 format UH0,29X, IPE12,5,/|X, iPtl2.5,5X, lPEl2»b,9X, 1PE12,53 
GO TO 2l6 
2/|/J LliSt=:l 

IF (IR1TE,EG, 1 3 GO TO 2 A' 2 


GQ ro 300 

216 CUMUvUE 
300 H = l^ AV C2)-?WAV(U 
Xlfvl50,0 
KNsN«2 

• DU 15 1 = 1, KN, 2 

XIMa;<INT+(iVir0)A('r.XU)+<),0ArtX(UU+wXU + 2}) 

5 ■ COKTIDUE 

IFfivruUX.EQ.n GO TU I'-tO’ 

XIKT 

H>-(2t-Rn = XlDT 

1/12 f'URwAT (///IX, =1P£12,5,8H XICHOWS) - 

fK = bC2 = 0, 

GO 10 146 
1/(0 XIM=X1NT/PI 
/jIRP = XlNT 
ldlKT = 0.0 
00 7 I=1,KN,2 

Y 1 MI s Y J + (H/3 ♦ 0 ) A ( W (; I ) + 4 , ( If I ) + A ( I + 2) ) 

7 CDNUNUH 

YIM = >! If^T/pr • 

EFX = XIN'UtnGRP/YTM 

TEFXsXIiNTarIGRP/ZIN'T 

OC=2,S*ALOGlO(l/EFX) 

BC2 = 2.5*AC.CG1 0 ( 1/TETX3 

144 00 126 1 = 1, fN 

126 i^AX( n = AX(I)*lsAV(I) 

X X I IV T = 0 . 0 
DO 126 1=1, KN, 2 

XXirvT = XXIM + (h/3,03A CWl'‘X(l) + 4,0AV<}WX(ttnfWyiX(It2)) 
12G CONTINUE 

yxiNTsXXiNT/FI 
EFFTsXXIN’I/XJNT 
IF(hFU-X,Eg, 13 GO TO 146- 
EFFlspUFFFT 
KHITC (6,1303 EFFT 

130 FiJkPAT (//lX,2iPLFFECTiVe WAVEIENGIH = lPEi2,5) 

146 CUNllMUE 

IF CNf-’UlJX»eG,0) I-AVOCFsEFFT 

iF(nfuux,er,03 go to a 

Z N 0 1'’ = 0 • ■ 

C Al U SPC n<y ( TF V F ( J ) , EFF T , ZNOn , bWA V , ZNON , ZNON, ZNON 3 
fJtN aV^UW AV/F'l 

IF ( J , ) ) GO 70 4451 . . 

H n A L F s b A V 


-p» 

•PS. 



AHREF = XlRP*Of^aGA 

PUK'CH 222, AHKEF 

?2? MJRF'ATCSX^APrFEFENCE FLUX s -frrlPfcia^S) 

PLINXH 223,WAVDFF,eW2£RC 

>23 F0RK-4T(5X,*QEF INED WAVELENGTH = A,n2,e,* ZERO BANDWIDTH = AfFT2. 

■ 16)' 

PUNCH 221 

>21 F'ORPATUX, *NO, TEFP WaVLNTH 8DW0TH HWTHOF BH HAG COL 

ihag in’t.flux view 
coMiNue 
■ Z N 0 N * 0 » 

call SPCTRM (TENP(U) ,'AAVDEF,ZN0N,BWAV0F fZNON,ZNON»ZNON) 

BwavdFsBwAVDF/PI 
IF(BV'AV,£Q,0, 9) GO TO 
dhv-.av=xInt/bwav 
DLWVDF = XlFT/f?WAVDF 
\i]in CONTINUE 

MAC:F AC = B>^'AV/eWREF 
>MA(;FAC = ALOGtO(MAGFAC) 

^ PAGF AC = "2.5*»'AGFaC 
• OHPAG = '2,SAALGGlC‘(XIHP*ontGA/AHHeF ) 

AlPP = XlHPAi|,5L-0fc LT 5 

Punch 200, icQL, TfP'^ CU9 ^EFFT ^ DtLWAV,OLWVUPrHAr.KACf BRHAGrXlRP, AIRP 
>00 F-UKNA7 ( IX, I3,2X,6FR,3,2E10, 5) 

If ( J-U) 6/8,8 
). CUkTTMJE 

IF (NrL'JX,fcG, U, GO TO 688 
NFl.UX=l , 

GO TO 887 
)88 CUNTINUE 

1F'(NN^N'TRANG)98,97,97 
97 JF(MJ-NFUNCT)95,96,96 

96 CUnTInUE 
STOP 
End 

subroutine SPCTRP (TEKP,NAVf RESP/W/WX/'aXX/HINTGI 

KnSERT DESlHfC E'MI3SJ.VITY HERE AS EMS s XXXX 

E N S 1 (( 0 0 

PI = 3, 1^1592659 

If (i-AV,EQ,0f I GO TO 10 

A 1 = 1 3679/ (wAVatL^^P/I 0000 ,) 

IF (AI.GT.68.) go to lO 
B1=EXP(AI) 

N=(3,7^l832E^l6)/( C (WAV/10000 *)aa5) A (BI-1 , ) ) aEMS 

go 10 12 



W=(3,7ai8 32E«i6;*CXP(«AI J / ( ( ( WA V / I OO 0 0 , ) ^ a5 ) ) 

WXsHK 

w X >: = >■' X / P I 

WiNTG=(3,66«6£-j[?) A (TP-'!p«*/4) /PI *cM3 

Rf TuRfJ 

fchO 
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A Point-Source Detection Routine 






The last routine presented here is a program developed for an early 
sky survey. The purpose in reviewing it here is to illustrate both 
the breadth of processing which can be done in a single pass bf the 
data and also the complexity of the required softV'/are. The routine 
unpacks and de-commutates the data and checks for errors and gaps. 
Three background channels are processed, the running noise 
computed, and the d^ta plotted. .Within the basic detection loop, 
the data is tested for signal peaks, correlation peaks, and signal 
length, and the radiation hits are separated from the data. Estimates 
of the amplitude and bias level are made^ and the position of the 
position of the signal is found from the time of detection. 


Inspection of a sample portion of the preliminary detection list 
reveals some of the basic problems which the following merging 
routines will need to deal with. The most complex problem is 
that the correlation coefficient does not track well with some of the 
other measurements of a good signal. For example, source number 22 
has a good correlation coefficient, but the estimated amplitude is 
less than half the peak height, and the amplitude estimate peak is 
significantly shifted in time from the data peak and also from the 
peak correlation coefficient. The correlation coefficient is below 
a reasonable error gate, but the same is true for signal number 10, 
where except for a slightly low p, the signal is very good. This 
pattern persisted throughout the data, and a careful study revealed 
that the- poor p values were the result of an uncertainty in the 
detector bias. Since the sensor used bidirectional logarithmic 
amplifiers, the uncertainty in bias led to a possible error in 
de-compressing the amplifier functions which would. tend to warp the 
signal shape significantly and degrade the value of the correlation 
coefficient, and also warping the noise spectrum. 


The software presented here v/as not designed to minimize its use 
of computer resources which would probably have resulted in a 
separation of the multiple functions of this routine. Furthermore, 
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using a maximum sensitivi^ty test which allowed a error rate, the 
resulting data was not significantly compressed. Of course, making 
multiple measurements ofWhe signal quality and not immediately testing 
on them. This allowed manual inspection of the data quality and 
careful adjustment of the tests which followed providing a sound study 
basis for a larger detection scheme. 
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INI TIALIZE DATA AHOOO PRELIMINARY CALCULATIONS 
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279 7 , 0 




ORIGINAL PAGE 
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^ rpTclTT? 9~ +': 2 ,"“+37"cX62 575 6+01,-0,11 5T0T5" - + 3 2 , + 0 . 1-1T2 pT 29£ + Q 2,-C ,499 
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. 5 + - 11 * 


__OAT_A {''HA\H I) r 5 C-.5) / C , C , 3 , ? , 0 , G i 3 , 1 3 C_2j 10 3_, 1 0 - , 10 & , 1 06 , 1 Q 
113 3,2c'’*', ll':',ill,il'2;:il3;illVil5VilV,117;ild , j ,‘20ir202,'203T2DV, 20 ’ 

Zr 5 - 2 i 2 ■, 21 : , 21 1 , 1 C= , 2 13 , 21- , 21 5 , 216 , 21 7 ,2 1 , 3 , 301 , 30 2, 50 
X3'3 >»3Gr, 3:‘6;3C7', 3 0 3", 2C ?. 511 , 312 , 31 S', 31 1 1 31 5 , 3 lil , 3 17,'31 & , 

“TTATa f FRTniriJTl ) I = 1 , 3J /lOhr^LZHANW 2, 13 Hi 21 SIG,10HNAL TRACES/ 

0A.TA"( ICH-NM n r j 'I = i ST’-l > / 6 , 9 , 1 0 i i' ,1 2 13 i ^rf 5 73 o 7l7Tr3''Tl 972 GT2T^' 2 2", ~ 
123,2^ , 26, 27 ,23,29, 30,31 ,32,i6,j-4,33,3c,37,35, 39,l3,4l,A2,A3,4-,-t5, 
297, ^S, •• 9, ;Q, -1, 72, 63,3-, 6 6, 5o, 57, 53,59,60, oil 6 2, 63/ 

“ivYAif <yv(i),i =r,‘2m : 

FO'-'MAT (3F25.14) 

”PTTT^T~l 6 ^ 

ORIMT -5, (I, YY(I} ,1=1,90 

■7^(TTT/rTTfX7rr{r7r277r4") 

, PRINT 16 



fSTAZ=- 23.7G 


u = c . L - « 

TJt}9=2G0. 



lio Y = 

00 21 1 = 1, t-D 


SUHY=SUMY+YY( I ) 


f u _i.ur = .'< 

SCNtND = SCNTIM + SCNLEN- 



CALL FLTID3 ( °ROGI D , ZL CNG , YMA X , ZFAC T ) 



QMS U ) = 0 . 


I Lrt ^ 











x?:ym=q. 












FSTR?C=PHR«:C 


















Y N H2 ( J ' = YWiTJ) ^ 
YNf;i(J)=Y?M J) t 

OTA {J,<) rY.M (j) 
SU’*i=sU''!-*OTA{J,K) 




IF RZFLVL ( J) =x?AR ( J) 

Tni5c-GTfTr^sRTjrYSTT;-smji 


IF { ({STGSG(J) “3.0^HSN( J) > .GT, :, ) .OR. (. NCT,.*iARKE?v{ J) ) ) 30 TO of. 




RMS(J>=3,^30RT (MSN (J) ) 



DTAKAX(J)=AMAX1(0TA(J,K) , OT AM AY ( J) ) 


GO TO 101 
TFT' 

GO TO iOl 


KSPB=K-»28 





. J. ) . 

KRJ(J)=I-23 


IJ J = h b U 

CONTINUE 


KAJ(J )=I-23 




COMTTNUE • 

AV“VAL( J) =REFLVL( J ) 


;N0(J)=PNS(J) 








MRJCT=NRJCT+1 


HRITHT3T~r)~ FHRi:CTJ,K'ST7?'r<( J) ,KrfA'KTT) , J I m.'I AX ( Jf, 3 GM Ifll J ) 

i,PKTlM( J) , ;«VRVAL( J ) .ZNO ( J) , NRJCT,RRPCNT ( J ) , KAJ (J) , amaxj ( J) ,KRJ (J) i 


o 

2RTr3XniT)”,r'TTAX(J) 

F0RMAT{5I5,r'E12.5,2I5,Z{l5,E12.5),Ei2.5) 

- 

c 

- 10-3 


c 

c 

rr— r-scu-RC'. m-s~ 32 EN ••Fm?:ro~KTTrcH~Excc:rDb ih£ referenc- l-vcil i 

THREE SIGMA, calculate THE OEGINING TIME AND THE TIM:: CF P-AK, 


"0 — 

AMj OUTPUl — rm.‘~Sl'A'K' UA 1 A TO fHz RECO^xDiNG t-iLz. 



IF ((. NCT. MARKER <J) ) .And. (K.EQ. 39)) KA J ( J) =KA J ( J) - 3 9 
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.U'/-<1 
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5 K ,'/nr n . <1 KT^U S'Q'BK-BKHN (TbkF^^BKP.N (I3K) ’ 

(n KSTGdRK) .LT.O.) BKSI G ( I BK ) = 0 . 
’KniGTl3K) = SQRT'r3KSlGTi 3K) ) 

'ONTIHUt 

r""i; fi'cr ; N^Fifisr; ^ 

'r-''-TN2=. false. 



nEL= 


I Li - L , U ' . -J . J 

QHSa )=SSRT (DEL) 


' » iUo c; -nr.to, IM , LOHA'N t i J , L i J , UMii ti ) , i = 0 , ) , 

1(1) ,CMS(I) tl = 27,4it) , (CHAM (I) ,KHS (I) »QMS ( I ) , I = 6 ,6 3 ) 



^4X,It-, ux»Ft . 2,ld ( 3 ( 


,13,-X,£10.-f,5X,E10.-)i/,19X 



1004 FORMAT (IHG, 5X,^TIME MAXIMUM REACHEo ’aT RECORD ^,I10.5X,+ FINAL TI 



c N ’ ) r i L - - 
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; 1 ■ ___/ 
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- 
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. IF(I0F(4)) S?,81 


= 
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GO TO ao 
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K-’-ii'.U / • 

.(“♦'1012J PHREC,TA, (CH^H(I) sRMSCI) ,OMS(I) ,1=8,23) ,(CHAH(I)' 

,RHS 

: rmTOT^STTTTT^ . , C 4 J , (IThTi N (T J ,hrs^m , 0Mb 1 1 J i 1 = - o i c ^ J 

IF(HGF(7)) 33,84 

ok 

‘A'- 1 Ti; { - , 1 ■: .l-ZT—gHRE C4'~A'rrUP7rrmTP^^TIT^, OMb U ) , I-O , 2p j , ioma rt i 

1(1) ,0^3(1) ,1=27,44), (CHAN(I) ,F^'S{I) ,QHS (I ) I = -6 ,53 ) 

t re I'l 

'.• 85 

G 0 T u ^ 

ENOFILE 4 ■ i 


T-C — 

• . • C ' 

NOW C0Ur4T UP THE OBJECTS D-TECTED, AND THE REJECTS, AND LIST T 

aSlES 


■r~C — = 

DO 2CC I=lt65 

















oecocG=oe:c^ 18 o , /pi 



































' Trr'rK — .rl 


4.0 APPENDIX; A SIGNAL PROCESSING GAME 


. The aim of this game' is to develop skills in signal processing. The. 

input data for this game are. the recorded data Uj^(t). It is assumed 

that the non-uniform scanning velocity has been corrected for already. 

The time coordinate is given in discrete numbers t = 0, 1, 2..;63. 

We may consider Uj(t) as being about one quarter of a single horizontal 

scan (e = constants). 

* ^ 

The rules of the games are as follov/s. The "investigator" gets the • 
sheet "Recorded Data Ug(t)" and the sheet "Problem #1." After 
solving this problem he Will give the solution to the "monitor" and' 
to the "game constructor." Now he may start- on problem #2, and so on. 
But it is' important that the investigator does not get the next 
problem sheet before he has finished the previous problem. The 
reason is that the formulation of the later problems contains parts- ■ 
of the answers to the earlier problems. This has to do with the 
basic structure of this simulation game: for performing any meaningful 

signal processing operation one must have some knov/ledge about the 
original signal and/or the noise. For example, in problem #1 the 
investigator is told that the noise is additive and non-negative. In 
the later problems, the investigator will be supplied with even 
more a priori information. Naturally, this should enable him to 
extract the signals better and better. But the methods for doing this 
increase in complexity. 

On the very last pages, following the problems, the design of the 
"recorded data" is explained, -and the true original signal is- 
unveiled. Obviously those pages should not be given to the investigator 
before he has solved all the problems. 
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. .RECORDED DATA Uj^(t) 


t is the discrete time variable running from t=0 to t=63’ 


t 

u,(t) 

t 

UrU) 

t 


t 

UR(t) 

0 
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16 
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48 

383 
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17 
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49 
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19 
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15 

51 

58 
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20 

20 

42 

36 
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52 

52 

5 
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09 

53 

66 

6 

07 

'22 

09 

38 

81 

54 

79 

7 

00 
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15 

39 

21 

55 

134 ■ 

'8 

62 

24 

13 
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81 

56 
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25 

■ 50 

41 

399 

57 
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10 
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26 

99 

42 

312 

58 
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11 
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27 

54 

.43 

348 

59 

94 

12. 

125 

28 

99 

■ 44 

303 

‘60 

56 

13 

125 

29 

35 

45 

383 

61 

67 

14 

173 

30 ' 

98 

■ 46 

317 

62 

51 

15 

00 

31 

02' 

47 

317 

63 
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Problem #1 


Given are the recorded data U|^(t) vrith t - 6^ 1...63. Wanted are the 
original data UQ(t), v/hich represent the “one-dimensional equivalent object 
radiation" S^(a). We assume that the known influences of the telescope 
[M(x',y'); R(x', y')] and of the electrical system [G(t)] have been 
compensated already or are negligible. But the recorded signal Uj^(t) is 
corrupted by additive noise M(t): 

- U^{t) = UgCt) + N(t) 

The only features known about the original signal UQ(t) and about the 
noise N(t) are that they are non-negative: 

UgCt) >. 0; N(t) >. 0. 

Furthermore, the noise N(t) is stationary, which means that the noise 
properties are not "drifting." In other words, short-term average 
•features of the noise remain the same from' the beginning to the end of 
the observation. 

Try to utilize the given a priori information for computing a new 
signal U-j(t) from U|^(t), v/hich somehow is better than Up(t) as an approximated 
representation of UgCt). Plot Ui(t) as a continuous curve, and also 
Uj^(t} .for comparison. 
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Problem #2 


Given are the facts: 

UQ{t)>0; N(t}>0; N = 50. 

By N we mean the linear average of the noise. This N can be visualized 
as the dark current of the photo receiver as measured with an instrument 
which rejects high frequencies. 

Based on these facts, try to compute a better signal UgCt) from 
U[^(t), Plot both U 2 (t) and U|^(t). 
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Problem #3 

I . 

Given ar^e the same facts as in the previous -problem. In addition, 
it is known that the noise- is approximately "whiite." 

' I 

N(t) = F + n(t); n{v) = |n(t) dv; 

V = m/64; m = -32, -31, ... -1, 0, +1., ... +30, +31; 

}fi(v)p constant. The amount of this "constant" is not known. Try 
to deduce it from the recorded data Ujj{t). You might have to make an 
intelligent guess. 
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Proble>/n #4 


Given are the same facts as in the previous problems j including the 
"constWt" which describes the noise power level. ^ 

|n(v)j2^i^^;rl0^ in - i < V < -f i 

Now that In(v)12 is known and Uj^( ) is computable, can you apply the- 
Wiener-filter theory, at least in a guessed approximation? Try it and 
compute U^(t). Plot U^(t) and Uj^(t). Hing: represent 1 Uq(v)I^ by 

a gaussian function of suitable peak power and width. Signal processing 
specialists always try it with a gaussian function if they don't know 
a better way. 



Problem §S 

Try the same approach as in the previous problem, but 
2 

guessed sine -shaped |Uq(v)|^ 

sinc^{v/vg); sine z = 

Plot the result Ug(t) and also Upj(t) for comparison. 


wi tH^ a 


sin^rz 

irz 
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Problem #6 


Try the same approach as in the previous problem, but a somewhat 
different guess for 1 Uq(v)|^ 

~ ^6 (Pq “ Pg) 

/ .1 

Herein 6 q means a function which isequal to 1 for v'- 0 and equal to 
0 for V / 0. Plot Ug(t) and Uj^(t). 


69 



Problem #7 

Based on all of the accumulated experience, try your own signal 
processing approach or simply guess what UgCt) might, have been. Call it 
U 7 (t). Plot UyCt) and Uj^(t). 
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